/*
 * Reaction.cpp
 *
 *  Created on: 23 May 2011
 *      Author: allan
 */

#include "Reaction.h"

// C++ includes
#include <cmath>
#include <iomanip>

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "NumericalUtilities.h"

Reaction::Reaction() :
kinetics(new ZeroKinetics())
{}

void Reaction::SetEquation(const string& reactant, const vector<string>& products, const vector<double>& stoichiometries)
{
	// Set the reactant species of the reaction
	this->reactant = {reactant, 0, 1.0};
	
	// Set the product species of the reaction
	for(unsigned i = 0; i < products.size(); ++i)
	{
		// Create a new product species
		Species product = {products[i], 0, stoichiometries[i]}; this->products.push_back(product);
	}
}

void Reaction::SetEquilibriumConstantData(const EquilibriumConstantData& eqConstant)
{
	this->eqConstant = eqConstant;
}

void Reaction::SetReactionKinetics(const ReactionKineticsPtr& kinetics)
{
	this->kinetics = kinetics;
}

void Reaction::SetIndexes(const SpeciesMap& speciesMap)
{
	// Set the index of the reactant species
	reactant.index = speciesMap[reactant.name];
	
	// Set the index of each product species
	BOOST_FOREACH(Species& product, products)
		product.index = speciesMap[product.name];
	
	// Set the index data of the kinetics object
	kinetics->SetIndexes(speciesMap);
}

const Reaction::Species& Reaction::GetReactant() const
{
	return reactant;
}

const vector<Reaction::Species>& Reaction::GetProducts() const
{
	return products;
}

const EquilibriumConstantData& Reaction::GetEquilibriumConstantData() const
{
	return eqConstant;
}

const double Reaction::EquilibriumConstant(double T, double P) const
{
	return eqConstant(T, P);
}

const double Reaction::ReactionQuotient(const VectorXd& a) const
{
	double Qr = 1.0/a[reactant.index];
	
	for(uint i = 0; i < products.size(); ++i)
	{
		const double ai = a[products[i].index];
		const double vi = products[i].stoichiometry;
		
		Qr *= std::pow(ai, vi);
	}
	
	return Qr;
}

const double Reaction::MassRate(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	return kinetics->MassRate(T, P, n, a);
}

ostream& operator<<(ostream& out, const Reaction& reaction)
{
	out << setw(10) << reaction.GetReactant().name << " (" << reaction.GetReactant().index << ")" << " = ";
	
	BOOST_FOREACH(const Reaction::Species& product, reaction.GetProducts())
		out << setprecision(2) << setw(4) << product.stoichiometry << setw(10) << product.name << " (" << product.index << ")" << ", ";
	
	return out;
}
